####Model result (table X10, model A29)####
main_cvm1t_t1_sb <- glm(as.formula(paste("cv_event_u", paste(c("sa_territory_t * adm_grp_rel_fragmentation","sa_territory_t * other_mean", unit_vars, nat_vars, "as.factor(region)","as.factor(year)","cv_event_u_peaceyears_l1","I(cv_event_u_peaceyears_l1^2)","I(cv_event_u_peaceyears_l1^3)","cv_event_slag_u_l1"), collapse = " + "), sep = " ~ ")), family=binomial(link="logit"),control = list(maxit = 50), data=main_unit)
main_cvm1t_t1_sb_cse <- data.frame(cluster.se(main_cvm1t_t1_sb, as.factor(main_unit$cowcode))[,2])
stargazer(main_cvm1t_t1_sb, se=c(main_cvm1t_t1_sb_cse), dep.var.labels.include = T, omit=c("cowcode|region|factor|years"), dep.var.labels = c("Non-state violence (unit)"), model.numbers =F, column.labels = c("(model A29)"), type="text", order=vars.order_t, title="Table X10. Additional results: unit-level analyses.", style = "apsr", star.cutoffs = c(.1, .05, .01, .001), star.char = c("†", "*", "**","***"), notes = c("† p<0.1; * p<0.05; ** p<0.01; *** p<0.001; country-clustered SE’s in parentheses.","The dependent variable is a binary variable equal to one if there is at least one instance of non-state violence in a given unit.","Region- and year-fixed effects included but not reported."), add.lines = list(c("Units included","all","min.-maj.","all","min.-maj.","all","min.-maj.")), notes.append=FALSE, covariate.labels = c("Territorial autonomy","Territorial autonomy x heterogeneity in unit","Territorial autonomy x size non-EPR groups", "Population (unit, logged)", "Area (unit, logged)", "Avg. ruggedness (unit)", "Oil in unit", "Distance capital (unit, logged)", "Distance border (unit, logged)","GDP pc. (logged)", "Population (state, logged)","Democracy","Ethnic fractionalization","Election year","Heterogeneity in unit","Size non-EPR groups","Civil violence spatial lag","Communal violence spatial lag"), out="../tables/additional_results_tables_appendices/tablex10.txt")
stargazer(main_cvm1t_t1_sb, se=c(main_cvm1t_t1_sb_cse), dep.var.labels.include = T, omit=c("cowcode|region|factor|years"), dep.var.labels = c("Non-state violence (unit)"), model.numbers =F, column.labels = c("(model A29)"), type="html", order=vars.order_t, title="Table X10. Additional results: unit-level analyses.", style = "apsr", star.cutoffs = c(.1, .05, .01, .001), star.char = c("†", "*", "**","***"), notes = c("† p<0.1; * p<0.05; ** p<0.01; *** p<0.001; country-clustered SE’s in parentheses.","The dependent variable is a binary variable equal to one if there is at least one instance of non-state violence in a given unit.","Region- and year-fixed effects included but not reported."), add.lines = list(c("Units included","all","min.-maj.","all","min.-maj.","all","min.-maj.")), notes.append=FALSE, covariate.labels = c("Territorial autonomy","Territorial autonomy x heterogeneity in unit","Territorial autonomy x size non-EPR groups", "Population (unit, logged)", "Area (unit, logged)", "Avg. ruggedness (unit)", "Oil in unit", "Distance capital (unit, logged)", "Distance border (unit, logged)","GDP pc. (logged)", "Population (state, logged)","Democracy","Ethnic fractionalization","Election year","Heterogeneity in unit","Size non-EPR groups","Civil violence spatial lag","Communal violence spatial lag"), out="../tables/additional_results_tables_appendices/tablex10.html")

####Figure A9 (appendix 3.6, non-mobilized groups)####
unit_cv_het <- margins_summary(main_cvm1t_t1_sb, variables = "sa_territory_t", at = list(adm_grp_rel_fragmentation = seq(0,1, length.out=25)), vcov = cluster.vcov(main_cvm1t_t1_sb, as.integer(subset(main_unit)$cowcode)))
unit_cv_other <- margins_summary(main_cvm1t_t1_sb, variables = "sa_territory_t", at = list(other_mean = seq(0,1, length.out=25)), vcov = cluster.vcov(main_cvm1t_t1_sb, as.integer(subset(main_unit)$cowcode)))
unit_cv_het_plot <- ggplot(data = unit_cv_het) + geom_line(aes(x=adm_grp_rel_fragmentation, y=AME)) + 
  geom_ribbon(aes(x=adm_grp_rel_fragmentation, ymin = AME - 1.96* SE, ymax = AME + 1.96*  SE), alpha = 0.2) +
  geom_hline(yintercept = 0, linetype = "dotted") + theme(plot.title = element_text(size=10)) +
  theme_bw() + theme(plot.title.position = "plot",legend.position = "bottom", text=element_text(family="Times"), axis.title=element_text(size=10), legend.text=element_text(size=10)) + 
  xlab('heterogeneity of unit') + ylab('change in predicted prob.\nof communal violence') + ggtitle("a) non-state violence, depending\non heterogeneity of unit") +
  scale_y_continuous(labels=scales::percent_format())
unit_cv_other_plot <- ggplot(data = unit_cv_other) + geom_line(aes(x=other_mean, y=AME)) + 
  geom_ribbon(aes(x=other_mean, ymin = AME - 1.96* SE, ymax = AME + 1.96*  SE), alpha = 0.2) +
  geom_hline(yintercept = 0, linetype = "dotted") + theme(plot.title = element_text(size=10)) +
  theme_bw() + theme(plot.title.position = "plot",legend.position = "bottom", text=element_text(family="Times"), axis.title=element_text(size=10), legend.text=element_text(size=10)) + 
  xlab("share of 'Others' in unit") + ylab('change in predicted prob.\nof communal violence') + ggtitle("b) non-state violence, depending\non share of 'Others' in unit") +
  scale_y_continuous(labels=scales::percent_format())
figurea9 <- grid.arrange(unit_cv_het_plot,unit_cv_other_plot, nrow=1, ncol=2)
ggsave(figurea9, file='../figures/figurea9.pdf', width = 18, height = 6, units="cm",dpi=1000)